library(funFEM)
library(reshape2)
library(ggplot2)
library(gridExtra)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(mclust)
library(stringr)
library(cluster)
library(clusterSim)
library(ggmap)
Cet ensemble de données contient des données provenant du système de partage de vélos de Paris, appelé Vélib. Les données sont des profils de chargement des stations de vélos sur une semaine. Les données ont été collectées toutes les heures pendant la période du dimanche 1er septembre au dimanche 7 septembre 2014.
On charge ici les données Velib disponibles dans la
librairie funFEM.
library(funFEM)
data(velib)
Ces données contiennent
data : les profils de chargement (nb de vélos disponibles / nb de points d’attache) des 1189 stations à 181 points de temps.
position : la longitude et la latitude des 1189 stations de vélos.
dates : les dates de téléchargement.
bonus : indique si la station est sur une colline (bonus = 1).
names : les noms des stations.
Afin d’avoir des semaines complètes, nous allons supprimer les 13 premières colonnes.
Velibdata <- velib$data[, -c(1:13)]
colnames(Velibdata) <- velib$dates[-c(1:13)]
On controle le nombre de jours et d’heures dans le jeu de données
day <- str_sub(colnames(Velibdata), 1, 3)
day <- factor(day, levels = c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim"))
table(day)
day
Lun Mar Mer Jeu Ven Sam Dim
24 24 24 24 24 24 24
hour <- str_sub(colnames(Velibdata), 5, 6)
table(hour)
hour
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
timeTick <- 24 * (0:6)
On stocke les informations sur les stations dans le dataFrame suivant :
station <- data.frame(nom = velib$names, velib$position, colline = velib$bonus)
On ajoute des informations complémentaires sur les stations :
# Les gares
gare <- rep(0, nrow(station))
gare[which(str_detect(station$nom, "GARE") == T)] <- 1
# Lieux touristiques (autour de Tour eiffel, Panthéon, Arc de Triomphe, Champs
# Elysées, Louvre, Notre Dame, Invalides, Trocadéro, Sacré Coeur, Montmartre)
tourisme <- rep(0, nrow(station))
TI <- c(209, 504, 897, 102, 518, 357, 256, 301, 992, 627, 579, 241, 646, 13, 283,
941, 976, 1189, 22, 633, 853, 973, 762, 195, 71, 131, 912, 302, 303, 1068, 774,
384, 307, 326, 436, 789, 250, 755, 461, 439, 918, 1065, 99, 151, 313, 871)
tourisme[TI] <- 1
# Les mairies
mairie <- rep(0, nrow(station))
mairie[which(str_detect(station$nom, "MAIRIE") == T)] <- 1
mairie[which(str_detect(station$nom, "HOTEL") == T)] <- 1
# Bois / Parc / Square
parc <- rep(0, nrow(station))
parc[which(str_detect(station$nom, "BOIS DE") == T)] <- 1
parc[which(str_detect(station$nom, "PARC") == T)] <- 1
parc[which(str_detect(station$nom, "SQUARE") == T)] <- 1
parc[which(str_detect(station$nom, "CANAL") == T)] <- 1
# Porte
porte <- rep(0, nrow(station))
porte[which(str_detect(station$nom, "PORTE") == T)] <- 1
station <- data.frame(station, gare = gare, tourisme = tourisme, mairie = mairie,
parc = parc, porte = porte)
Fonction pour tracer les profils (charge de chaque station / loading) en fonction des jours et heures pour les stations choisies.
plotprofils <- function(Velibdata, station, numstation, plotsem = TRUE) {
# Velibdata = les données de velib initiales station = ens. des
# informations sur les stations numstation = vecteur contenant le numéro de
# ligne des stations dont on souhaite tracer les profils
library(reshape2)
library(ggplot2)
Dataaux <- data.frame(id.s = station$nom, Velibdata)
Aux <- melt(Dataaux[numstation, ], id = c("id.s"))
Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
5, 6)))
Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
"Ven"))
if (plotsem == TRUE) {
g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
ncol = 2) + xlab("Time") + ylab("Loading") + scale_x_continuous(breaks = seq(0,
24, 4)) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
labs(col = "Classif.") + scale_x_continuous(breaks = c(0, 7, 9, 12, 14,
20, 24))
} else {
v <- rep(0, nrow(Aux))
for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
levels(Aux$day)[j])] + (24 * (j - 1))
Aux <- data.frame(Aux, v = v)
rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
ymax = rep(1, 7), color = rainbow(n = 7))
g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
ylab("Loading") + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
}
return(g)
}
Exemple d’utilisation :
plotprofils(Velibdata, station, numstation = c(813, 655, 468))
plotprofils(Velibdata, station, numstation = c(813, 655, 468), plotsem = FALSE)
Se donnant une classification, fonction pour tracer les profils moyens par classe :
plotprofilsmoy <- function(Velibdata, clustering, plotsem = TRUE) {
library(reshape2)
library(ggplot2)
Dataaux <- matrix(0, nrow = max(clustering), ncol = ncol(Velibdata))
for (k in 1:max(clustering)) {
Dataaux[k, ] <- apply(Velibdata[which(clustering == k), ], 2, mean)
}
colnames(Dataaux) <- colnames(Velibdata)
Dataaux <- data.frame(id.s = as.factor(1:max(clustering)), Dataaux)
Aux <- melt(Dataaux, id = c("id.s"))
Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
5, 6)))
Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
"Ven"))
if (plotsem == TRUE) {
g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
ncol = 2) + xlab("Time") + ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") +
theme(axis.title.x = element_blank()) + labs(col = "Classif.") + scale_x_continuous(breaks = c(0,
7, 9, 12, 14, 20, 24))
} else {
v <- rep(0, nrow(Aux))
for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
levels(Aux$day)[j])] + (24 * (j - 1))
Aux <- data.frame(Aux, v = v)
rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
ymax = rep(1, 7), color = rainbow(n = 7))
g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
}
return(g)
}
Exemple d’utilisation avec les stations sur colline par rapport aux autres :
plotprofilsmoy(Velibdata, clustering = station$colline + 1, plotsem = TRUE)
plotprofilsmoy(Velibdata, clustering = station$colline + 1, plotsem = FALSE)
Fonction auxiliaire pour comparer une classification avec les informations auxiliaires disponibles pour les stations
stationcaract <- function(station, clustering) {
library(dplyr)
# Dataaux<-data.frame(id.s=c(1:nrow(Data)),station) aux <- cbind(Dataaux,
# clust)
aux <- cbind(station, clustering)
aux.long <- melt(data.frame(lapply(aux[, -c(2:3)], as.character)), stringsAsFactors = FALSE,
id = c("nom", "clustering"), factorsAsStrings = T)
# Effectifs
aux.long.q <- aux.long %>%
group_by(clustering, variable, value) %>%
mutate(count = n_distinct(nom)) %>%
distinct(clustering, variable, value, count)
# avec fréquences
aux.long.p <- aux.long.q %>%
group_by(clustering, variable) %>%
mutate(perc = count/sum(count)) %>%
arrange(clustering)
gaux <- ggplot(data = aux.long.p, aes(x = clustering, y = perc, fill = value)) +
geom_bar(stat = "identity") + facet_wrap(~variable)
return(gaux)
}
On reprend ici rapidement quelques éléments d’analyse descriptive vus précédemment dans l’UF.
library(corrplot)
timeTickAux <- c(timeTick, 168)
for (k in 1:7) corrplot(cor(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k + 1])]),
method = "ellipse")
velibmeanday <- matrix(0, nrow = nrow(Velibdata), ncol = 7)
for (k in 1:7) {
velibmeanday[, k] <- apply(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k +
1])], 1, mean)
}
colnames(velibmeanday) <- c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim")
ggplot(melt(velibmeanday), aes(x = Var2, y = value)) + geom_boxplot() + xlab("") +
ylab("Loading")
velibHour <- data.frame(value = colMeans(Velibdata), day = day, hour = as.numeric(hour))
ggplot(velibHour, aes(x = hour, y = value, col = day)) + geom_line() + geom_point() +
ylab("Loading") + ggtitle("Hourly loading, averaged over all stations")
library(ggmap)
# Average loading for each station
load <- rowMeans(Velibdata)
qmplot(longitude, latitude, data = station, color = load) + scale_colour_gradient(high = "red",
low = "yellow")
# scale_colour_gradient(high='blue',low='white')
# Loading moyen à l'heure t (le premier temps est 0h)
t <- 0
loadheure <- rowMeans(Velibdata[, seq(from = (t + 1), by = 24, length = 7)])
qmplot(longitude, latitude, data = station, color = loadheure) + scale_colour_gradient(high = "red",
low = "yellow")
qmplot(longitude, latitude, data = station, color = as.factor(station$colline))
qmplot(longitude, latitude, data = station, color = as.factor(station$tourisme))
Question : Faites une ACP des stations et analysez les résultats.
respca <- PCA(.......)
fviz_eig(....)
fviz_pca_var(......)
fviz_pca_ind(........)
Correction :
library(FactoMineR)
library(factoextra)
respca <- PCA(Velibdata, scale.unit = T, graph = F, ncp = 15)
# Variance expliquée
fviz_eig(respca)
ggplot(data.frame(comp = 1:20, cumul = respca$eig[1:20, 3]), aes(x = comp, y = cumul)) +
geom_line() + geom_point()
# Visualisation des variables
fviz_pca_var(respca, axes = c(1, 2), geom = c("arrow"))
fviz_pca_var(respca, axes = c(2, 3), geom = c("arrow"))
for (k in 1:7) corrplot(respca$var$cor[(timeTickAux[k] + 1):(timeTickAux[k + 1]),
1:3], method = "ellipse")
# Projection des individus
fviz_pca_ind(respca, axes = c(1, 2), geom = c("point"), habillage = as.factor(station$colline))
fviz_pca_ind(respca, axes = c(1, 3), geom = c("point"))
nbdim <- 4
aux <- melt(respca$var$coord[, 1:nbdim])
aux$Var1 <- str_sub(aux$Var1, 1, 3)
aux <- data.frame(aux, time = rep(c(0:167), nbdim))
rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(min(aux$value),
7), ymax = rep(max(aux$value), 7), color = rainbow(n = 7))
ggplot(data = aux) + geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin,
ymax = ymax, fill = color), alpha = 0.3) + geom_line(data = aux, aes(x = time,
y = value)) + facet_wrap(~Var2, ncol = 2) + geom_hline(yintercept = 0, col = "red") +
xlab("") + ylab("") + scale_x_continuous(breaks = seq(0, 168, 24)) + theme(legend.position = "none")
Dans cette partie, nous allons utiliser les coordonnées de l’ACP comme matrice de données pour classer les stations.
velibacp <- respca$ind$coord[, 1:7]
Question : A l’aide d’une procédure des Kmeans, déterminez une classification des stations.
# A COMPLETER
Correction :
set.seed(12345)
Kmax <- 20
reskmeans <- matrix(0, nrow = nrow(velibacp), ncol = (Kmax - 1))
Iintra <- NULL
Silhou <- NULL
for (k in 2:Kmax) {
resaux <- kmeans(velibacp, k, nstart = 20, iter.max = 30)
reskmeans[, (k - 1)] <- resaux$cluster
Iintra <- c(Iintra, resaux$tot.withinss)
aux <- silhouette(resaux$cluster, daisy(velibacp))
Silhou <- c(Silhou, mean(aux[, 3]))
}
rm(resaux, aux)
df <- data.frame(K = 2:Kmax, Iintra = Iintra, Silhou = Silhou)
g1 <- ggplot(df, aes(x = K, y = Iintra)) + geom_line() + geom_point() + xlab("Nombre de classes") +
ylab("Inertie intraclasse")
g2 <- ggplot(df, aes(x = K, y = Silhou)) + geom_line() + geom_point() + xlab("Nombre de classes") +
ylab("Critère Silhouette")
grid.arrange(g1, g2, ncol = 2)
resICL <- mclustICL(velibacp, G = 1:20, modelNames = c("EII"))
summary(resICL)
Best ICL values:
EII,11 EII,13 EII,19
ICL -41830.01 -41855.1410 -41862.86314
ICL diff 0.00 -25.1267 -32.84888
Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue.
Correction :
On retient la classification en \(K=5\) classes :
Kkmeans <- 5
classifkmeans <- reskmeans[, Kkmeans - 1]
table(classifkmeans)
classifkmeans
1 2 3 4 5
406 151 209 241 182
On visualise la classification sur la carte 2D :
qmplot(longitude, latitude, data = station, colour = as.factor(classifkmeans))
Profil moyen par classe :
plotprofilsmoy(Velibdata, clustering = classifkmeans, plotsem = TRUE)
plotprofilsmoy(Velibdata, clustering = classifkmeans, plotsem = FALSE)
Croisement avec les informations auxilaires sur les stations :
stationcaract(station, clustering = classifkmeans)
Visualisation d’une classe par exemple
I <- which(classifkmeans == 5)
qmplot(longitude, latitude, data = station[I, ])
Question : A l’aide d’une procédure hiérarchique, déterminez une classification des stations.
# A COMPLETER
Correction :
Kmax <- 20
resward <- hclust(dist(velibacp, method = "euclidean"), method = "ward.D2")
df <- data.frame(K = 1:Kmax, height = sort(resward$height, decreasing = T)[1:Kmax])
ggplot(df, aes(x = K, y = height)) + geom_line() + geom_point()
CH <- NULL
Kmax <- 20
for (k in 2:Kmax) {
CH <- c(CH, index.G1(velibacp, cutree(resward, k)))
}
daux <- data.frame(NbClust = 2:Kmax, CH = CH)
ggplot(daux, aes(x = NbClust, y = CH)) + geom_line() + geom_point()
Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par CAH. Comparez cette classification avec celle obtenue par les Kmeans.
# A COMPLETER
Correction :
Les critères précédents proposent de retenir 3 ou 5 classes :
classifCAH3 <- cutree(resward, 3)
classifCAH5 <- cutree(resward, 5)
table(classifCAH3, classifCAH5)
classifCAH5
classifCAH3 1 2 3 4 5
1 163 0 0 0 180
2 0 237 180 0 0
3 0 0 0 429 0
On visualise la classification sur la carte 2D :
qmplot(longitude, latitude, data = station, colour = as.factor(classifCAH5))
Profil moyen par classe :
plotprofilsmoy(Velibdata, clustering = classifCAH5, plotsem = TRUE)
plotprofilsmoy(Velibdata, clustering = classifCAH5, plotsem = FALSE)
Croisement avec les informations auxilaires sur les stations :
stationcaract(station, clustering = classifCAH5)
Visualisation d’une classe par exemple
I <- which(classifCAH5 == 2)
qmplot(longitude, latitude, data = station[I, ])
Comparaison avec la classification des Kmeans :
table(classifkmeans, classifCAH5)
classifCAH5
classifkmeans 1 2 3 4 5
1 0 22 3 381 0
2 0 119 26 0 6
3 12 37 0 0 160
4 151 33 0 43 14
5 0 26 151 5 0
adjustedRandIndex(classifkmeans, classifCAH5)
[1] 0.6411121
Question : A l’aide des mélanges gaussiens, déterminez une classification des stations.
# A COMPLETER
Correction :
set.seed(12345)
resICLall <- mclustICL(velibacp, G = 2:20)
summary(resICLall)
Best ICL values:
VEV,7 VVE,10 VVE,7
ICL -40227.06 -40271.08062 -40281.27660
ICL diff 0.00 -44.02224 -54.21822
resICL <- Mclust(velibacp, G = 7, modelNames = "VEV")
resBICall <- mclustBIC(velibacp, G = 2:20)
summary(resBICall)
Best BIC values:
VVE,10 VEI,13 VEE,13
BIC -39953 -39975.77232 -39981.70915
BIC diff 0 -22.77722 -28.71404
resBIC <- Mclust(velibacp, G = 10, modelNames = "VVE")
Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par modèles de mélange. Comparez cette classification avec celles obtenues précédemment.
# A COMPLETER
Correction :
Aux <- data.frame(label = paste("Cl", resICL$classification, sep = ""), proba = apply(resICL$z,
1, max))
h1 <- ggplot(Aux, aes(x = label, y = proba)) + geom_boxplot()
h2 <- fviz_cluster(resICL, data = velibacp, ellipse = F, geom = "point") + ggtitle("") +
theme(legend.position = "none")
grid.arrange(h1, h2, ncol = 2)
plotprofilsmoy(Velibdata, clustering = resICL$classification, plotsem = TRUE)
plotprofilsmoy(Velibdata, clustering = resICL$classification, plotsem = FALSE)
qmplot(longitude, latitude, data = station, colour = as.factor(resICL$classification))
table(resICL$classification)
1 2 3 4 5 6 7
331 186 83 207 132 128 122
Comparaison entre la classification avec ICL et celle des Kmeans
table(resICL$classification, classifkmeans)
classifkmeans
1 2 3 4 5
1 42 89 79 89 32
2 44 0 0 142 0
3 2 0 0 0 81
4 190 0 0 8 9
5 0 0 130 2 0
6 128 0 0 0 0
7 0 62 0 0 60
Dans cette partie, on va s’intéresser à un clustering sur un jeu de données “résumé” au sens suivant : on fait la moyenne des loadings selon les tranches horaires suivantes : 0-7h, 8h-20h et 21h-23h pour chaque jour et chaque station.
VelibMean <- matrix(0, nrow = nrow(Velibdata), ncol = 3 * 7)
for (jour in 1:7) {
VelibMean[, 3 * (jour - 1) + 1] <- apply(Velibdata[, 24 * (jour - 1) + c(1:8)],
1, mean)
VelibMean[, 3 * (jour - 1) + 2] <- apply(Velibdata[, 24 * (jour - 1) + c(9:21)],
1, mean)
VelibMean[, 3 * (jour - 1) + 3] <- apply(Velibdata[, 24 * (jour - 1) + c(22:24)],
1, mean)
}
colnames(VelibMean) <- paste(rep(unique(day), each = 3), "-", rep(c("0-7h", "8h-20h",
"21-23h"), 7), sep = "")
Question : Etudiez les corrélations entre variables
du jeu de données VelibMean.
# A COMPLETER
Correction :
corrplot(cor(VelibMean[, ((1:7) * 3) - rep(c(2, 1, 0), each = 7)]), method = "ellipse")
Question : Faites une ACP et commentez.
# A COMPLETER
Correction :
resacpMean <- PCA(VelibMean, scale.unit = T, graph = F, ncp = 10)
# Variance expliquée
fviz_eig(resacpMean)
# Visualisation des variables
fviz_pca_var(resacpMean, axes = c(1, 2), geom = c("arrow", "text"))
fviz_pca_var(resacpMean, axes = c(2, 3), geom = c("arrow", "text"))
corrplot(resacpMean$var$cor[, 1:3], method = "ellipse")
# Projection des individus
fviz_pca_ind(resacpMean, axes = c(1, 2), geom = c("point"), habillage = as.factor(station$colline))
fviz_pca_ind(resacpMean, axes = c(1, 3), geom = c("point"))
Question : Faites une classification des stations à
partir du jeu de données VelibMean à l’aide des Kmeans.
Etudiez la classification obtenue.
# A COMPLETER
Correction :
set.seed(12345)
Kmax <- 20
reskmeansbis <- matrix(0, nrow = nrow(VelibMean), ncol = (Kmax - 1))
Iintra <- NULL
for (k in 2:Kmax) {
aux <- kmeans(VelibMean, k)
reskmeansbis[, (k - 1)] <- aux$cluster
Iintra <- c(Iintra, aux$tot.withinss)
}
df <- data.frame(K = 2:Kmax, Iintra = Iintra)
ggplot(df, aes(x = K, y = Iintra)) + geom_line() + geom_point()
Etude de la classification obtenue
classifkmeansbis <- reskmeansbis[, 4]
table(classifkmeansbis)
classifkmeansbis
1 2 3 4 5
174 254 200 253 308
plotprofilsmoy(Velibdata, clustering = classifkmeansbis, plotsem = FALSE)
plotprofilsmoy(Velibdata, clustering = classifkmeansbis, plotsem = TRUE)
qmplot(longitude, latitude, data = station, colour = as.factor(classifkmeansbis))
Comparaison avec la classification obtenue sur les coordonnées de l’ACP :
table(classifkmeans, classifkmeansbis)
classifkmeansbis
classifkmeans 1 2 3 4 5
1 1 5 121 0 279
2 9 98 0 44 0
3 8 0 0 201 0
4 156 1 76 8 0
5 0 150 3 0 29
Comparaison avec les variables colline et tourisme
table(classifkmeansbis, station$colline)
classifkmeansbis 0 1
1 172 2
2 253 1
3 149 51
4 252 1
5 236 72
table(classifkmeansbis, station$tourisme)
classifkmeansbis 0 1
1 174 0
2 236 18
3 199 1
4 248 5
5 286 22
Question : Faites une classification des stations à
partir du jeu de données VelibMean à l’aide des mélanges
gaussiens. Etudiez la classification obtenue.
# A COMPLETER
Correction :
set.seed(12345)
resICLallbis <- mclustICL(VelibMean, G = 10:20, modelNames = c("VVE", "VEV", "EVV",
"VVV"))
summary(resICLallbis)
Best ICL values:
VVE,15 VVE,16 VVE,14
ICL 22181.01 22079.2176 22070.926
ICL diff 0.00 -101.7971 -110.089
resICLbis <- Mclust(VelibMean, G = 16, modelNames = "VVE")
Aux <- data.frame(label = paste("Cl", resICLbis$classification, sep = ""), proba = apply(resICLbis$z,
1, max))
h1 <- ggplot(Aux, aes(x = label, y = proba)) + geom_boxplot()
h2 <- fviz_cluster(resICLbis, data = velibacp, ellipse = F, geom = "point") + ggtitle("") +
theme(legend.position = "none")
grid.arrange(h1, h2, ncol = 2)
table(resICLbis$classification)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
74 192 92 83 53 125 54 63 52 34 43 54 37 105 54 74
plotprofilsmoy(Velibdata, clustering = resICLbis$classification, plotsem = FALSE)
plotprofilsmoy(Velibdata, clustering = resICLbis$classification, plotsem = TRUE)
qmplot(longitude, latitude, data = station, colour = as.factor(resICLbis$classification))
stationcaract(station, resICLbis$classification)
Les données de Vélib sont des données fonctionnelles (on étudie un
phénomène qui évolue au cours du temps). Il existe des méthodes de
clustering dédiées aux données fonctionnelles, comme par exemple le
package funFEM.
Pour utiliser cette méthode, on commence par considérer la décompostion dans la base de Fourier des courbes.
basis <- create.fourier.basis(c(0, 167), nbasis = 25)
fdobj <- smooth.basis(1:168, t(Velibdata), basis)$fd
La méthode funFEM est une procédure de clustering basée
sur des mélanges fonctionnels. Comme pour les mélanges gaussiens, il y a
différentes formes disponibles selon les contraintes imposés dans le
modèle (le détail est hors programme !). Cette méthode est détaillée
dans l’article de Bouveyron et al. (2015).
Question : Executez les commandes suivantes et exploitez les résultats. On considère une classification en 10 classes comme dans Bouveyron et al. (2015).
resfunFEM <- funFEM(fdobj, K = 10, model = "AkjBk", init = "kmeans", lambda = 0,
disp = TRUE)
>> K = 10
AkjBk : bic = -20459.07
The best model is AkjBk with K = 10 ( bic = -20459.07 )
df <- data.frame(proba = apply(resfunFEM$P, 1, max), classif = as.factor(apply(resfunFEM$P,
1, which.max)))
ggplot(df, aes(x = classif, y = proba)) + geom_boxplot()
classifFEM <- apply(resfunFEM$P, 1, which.max)
plotprofilsmoy(Velibdata, clustering = classifFEM, plotsem = FALSE)
plotprofilsmoy(Velibdata, clustering = classifFEM, plotsem = TRUE)
qmplot(longitude, latitude, data = station, colour = as.factor(classifFEM))
stationcaract(station, classifFEM)
table(resICL$classification, classifFEM)
classifFEM
1 2 3 4 5 6 7 8 9 10
1 16 59 30 24 24 7 2 48 0 121
2 13 3 117 0 0 0 37 15 1 0
3 0 0 0 0 0 70 0 13 0 0
4 0 0 1 0 0 1 83 88 34 0
5 59 1 0 70 0 0 0 0 0 2
6 0 0 0 0 0 0 4 8 116 0
7 0 0 0 0 65 34 0 2 0 21
adjustedRandIndex(resICL$classification, classifFEM)
[1] 0.3506168
table(resICLbis$classification, classifFEM)
classifFEM
1 2 3 4 5 6 7 8 9 10
1 0 2 7 0 0 0 21 40 0 4
2 10 2 103 0 0 0 42 24 10 1
3 0 0 0 13 63 4 0 1 0 11
4 0 1 0 0 9 18 0 21 0 34
5 0 0 0 0 5 43 0 5 0 0
6 14 4 9 24 1 2 0 3 0 68
7 0 0 1 0 0 0 50 1 2 0
[ getOption("max.print") est atteint -- 9 lignes omises ]
adjustedRandIndex(resICLbis$classification, classifFEM)
[1] 0.2955738